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ome UN CRODUC TION 


It is frequently desired to generate .Poisson random 


variates in simulations. There are standard exact methods 


for doing this; the problem arises when a computer is used 


to generate the Poisson random number which has a large mean. 
For example, generating one random number such as 105 "from 2 
Poisson distribution with mean 100 needs at least 105 calls 
to a pseudo-random number (uniform (0,1)) generator. Compu- 
ter time requirements become important cost factors when 
considering various methods for generating random numbers. 

The objective of this paper is to examine several ap- 
proximated ways of generating Poisson random variables and 
to determine the method which gives minimum execution time 
and small mean absolute deviation according to the Poisson 
mean value. The mean value is the only parameter in the 
disitribution. Comparison statistics to determine the best 
approximation to the Poisson distribution are the cumulative 
probability, mean absolute deviation and Kolmogorov-smirnov 
test. Finally a composite generating procedure according to 
The mean value is suggested. 

The following notation is used in this study. 

U, denotes a uniform (0,1) random variable; 
N denotes a Poisson random variable where mean is m; 


Z denotes a random variable from a unit normal distribution. 





II. GENERATION OF POISSON DISTRIBUTED VARIATES 


i. een OLSSON DES TRIBU PION 
A random variable N with integer values has a Poisson 


suri buLLon if 





prob {N = -m} = am = @sis2Ponn « 


In order to generate a Poisson random number N from a Poisson 
distribution with mean m, the following algorithm is pre- 
sented. It is.the standard exact method for generating these 
deviates. 

Let Us i =1,2,..., be independent uniform (0,1) random 


Voariares meine otsson variave, N, is computed as: 


where N is non distributed as a Poisson with mean m, i.e., 
Prob (N=n) = e 'm™/n! 
Equivalently, since the logarithm is a monotone transforma- 


tion, we have 


Ogle 1n(U, ) < In(e ”) = «=m 





oS 

k k ae 
Cetin foe) = mac.) > 1h Come 

i“ i 76 ih 

i=1 1=1 
| k+1 k+1 
=z & In(U, ) = In ( Il U,). 
L=1 j=] 





Changing the signs and direction of the inequalities we 


have: 


0 if -1n(U,) >m 


N = 
k kt1 
k if E - In(U,)<m< £ - In(U,) 
i=1 Peal = 
or letting E. = = In(U, ) 
jie 18 > in 
L- 
i 
k kt1 
edie ees EB. <m < i Be 
i=l i=l 





where the EB, are exponentially distributed variates with 
mean aly 

If n multiplications of uniform (0,1) random numbers is 
strictly greater than e" and if ntl MULCL pla Ca tl Cneem sun 1 ior 
(0,1) random number is equal and less than e " then n is the 
Poisson random number. Generally generating one random num- 
ber from Poisson process with parameter m requires on the 
average mt+tl uniform (0,1) random numbers. This is because 
the number generated is ntl. When m is large it is clear 
that generating Poisson random numbers with the above method, 
although it is exact takes a lot of computer time and this 
Logioumidvwepe UNeCOMoml Cat. =! Add1GIOmM, thes large number 
of multiplications can produce serious precision problems 


om ag dagi val computer. 
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B. APPROXIMATIONS FOR POISSON VARIATES 
1. Normal Approximation | 

In a Poisson process with parameter i» it is necessary 
to Eenerate random variables from the Poisson distribution 
with parameter (mean) m. Now look at counts in (0,x), where 
x satisfies Ax =m. The central limit theorem says that as 
m goes to positive infinity, (or when x gives the infinity in 
Poisson process with fixed 4), then N, which has mean m and 
variance m, is such that Neen is approximately distributed 
ads a2 unit normal random variable. Denote a random variable 
Preonea Unutlenormaledistribution by Z. so N Lswoastriplwkeed 
approximately as m *z +m. In order to generate Poisson ran- 
dom numbers from the normal distribution, first generate 2Z; 
voldteig) 1hea2 


7 
0 SMA, ar ie, i 


2 
It 


15 
{ m Z++tm- ee | otherwise 


where fal denotes theygreatest integer less than or equal to 
a; also known as the "floor" of a. N is then the approximated 
Poisson random variate. 
2. square Root ‘Transformation of Poisson Distribution 
If N is a Poisson random variable with mean m then 
Y = YN + 378 is approximately distributed as a normal dis- 
“6 


tribution with mean m°* and variance *. This result is due to 


Bartlett [Ref. 9], 


aE 





This method is derived as follows: let Y = YN +C ~ N(u,07) 
where C is a non-negative constant. Let t = N =m and 
ieee nen OC. Define coefiicients for s = 1g2,3,...by 


@ = (-1)8tL Le(-1)+(-3)+++(-28+3) 
. 2° ss! 


Pichi etOre any el en Mie weeiWave=aglaylor series expansion. 


r 
S 1, 


KK 
tl 


iF t 2 S t 
t aime ——— eee -_ —_——_ 
vm {1+A, = A, (7+) + +(-1) Asaa 


+ Roe \ 


wal 
Ifet > 0. wels@e ,at once |R. | < Axta/ (a's 2 converges and 
is bounded (F. J. Anscombe (4)). We note now that the moments 


of t are uw, = 0, Wp =m, U3 =m, Wy = 3m* + m,..., which give 


SaG ., SOUS BAU ae ae 


Var (Y) ~ % (1 + oe 3amn2 ye 


so that when C = 3/8 Var (Y) ~ % (1 + 1/16m?2). Also 


Le 1 Bye 
E(Y) ym+C — —z— + =~ + 
8m = 128m’ 


~ 


Let XNR = YN + 3/8. Then XNR is approximately normally dis- 


tributed with mean vm and variance %. 


XNR - vm + 3/8 


- i 
XNR = £ + vn + 378 , 
WN ¥ 378 = S + ym + 378 , 
thus set 


al 





0 if (E+ Vmt37B) - 3/8 < l 


Ne 
(3 + /m+3/8) 7- 3/8 | . otherwise 
N is then the approximated Poisson variate. 


We now need to calculate the probability distribution 
of N obtained in this way from the square root transformation. 


We want the probability that 


n - 1 + 3/8 < XNR* <n + 3/8 
if we divide by the variance *% we get 
K(n - 1 + 3/8) < XNR? 4 < A(n + 3/8). 
Note that 4-XNR* is distributed with a non-central x{ distri- 


bution. The non-central x{ density is 


-(xtu?/o7) - 4 u 
f(x) = a eS ovX + ehot = ] 
ss 2v¥x Von 


1 
Thus if wsm-, o%=(%)2, then 


eM) (Mt) VE gg (YM/) 


f (x) = 
. 2¥x yom 
and 
N+3/8 N+3/8-1 
(45) ¢ (35) 2 
prob(x=N) = f(x) dx _ f(x) dx 


O 
Witsoe Nilowss Us8vo valvatel directly how well the 
G2 Sle PoUplOnwoheN approximates the distribution of a ‘Poisson 


variate with parameter m. 


ILS 





Note that since in the LLRANDOM package it takes the 
same Serie Cl surme co peneraves5 uniform random variables 
as it takes to generate a normal random var slabaed: the pro- 
cedures will be competitive timewise once m is much greater 
ean ee 

3. Cube Root Transformation of Del eson Drs CempuL 1 On 

If N is a Poisson random variable with mean m then. 

Y = YN-1/24 is approximately distributed as a normal distri- 


1/3 and variance 1/9°vm when N > 1. This 


bution with mean m 
comes from the following derivation which is essentially the 
same =proceduremas Tor the® squidre root transformation. Sup— 


3 
pose Y = YNtC is distributed as N(u,o7). Let t = N - m and 


m' =m +C. Define coefficients for s = 1,2,3,... by 
aw £2874 1 (-2) (-5) (-8) ++ 38) 
Ss 38 .i : 


ita) eee 1 We have The laylowsseries expansien 
3 
C i Gy 3 Cyy 
a= ’ 7 —— ew —e —_ — — ® ® 
Meanie (1 + ja (Sp) a3 (=7) ay (—+) 


oot (al earl (+ S-li aR 


s—1 S 


if t > 0, we see at once that IR. ,t,/¢m' 8073 converges. 


Therefore, 


a S T,s-—l 
i +(-1) ae 


»\-l/3 _ eel /3 
R, (m ) =(l1+— ) = (at &e a. a 


} 





R.(mt)~7/3 7 
—— = 


bounded. 


obey ee 


Ss 


converges and is 


ttm 8 


4 


We note now that the moments of t are u,=0, U,=m, 


Les Wiel THe eae Saving us 


ee 


E(Y) 


and 


He 


Var(Y) 


If C = =o then 


var(Y) 


Let 
YNR = 
a 
YNR = 
VY N=4ed = 
Thus set 
0 











wy aaa 
~¥m+tC - a3) "373 oo 


mM 


1 2m*+m ei 


IL ale 
Cov)? fo carye eM - BT tye 


1 2), ets cage 





it 
9 ° vin 





°YN - 1724 N> 1 


VN / nied. Ut 


mL: 
— :3uvn 





3/m + 1/24 +(1/3°vm) Z 


3 /m tom/2h + (1/3'vm) Z. 


if (3v¥mt1/2) 5 6 /m Z) or <call 


il. 
cae Vint l/ 24 + 3 ®Ym 7)? 4+ sir otherwise 





N is now the approximated Poisson variate. 
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III. EVALUATION OF THE METHODS 


Generally, generating Poisson random number from the 
exact method is known to take a long execution time, since 
one generated random number "N" requires on the average N 
multiplications of uniform (0,1) ot oe Hohn erase Therefore, 
epenerating Poisson variables with the exact method (which 
gives the best accuracy) is good for a small mean, m, while 
Hic. approx iIntclon methods, which Cake shorter Seyarteeys vie 
but with less accuracy, should be used for large m. ene we 
‘need a trade-off between execution time and accuracy to 
choose the generation method according to the mean value of 
Pic wae? S SoMa! sori t om. 

The comparison statistics show how closely the methods 
approximate the original Poisson distribution. In the 
Kolmogorov-Smirnov test, all approximations are accepted at 
Significance level a = .05. From the comparison of the em- 
Ditrical probabilicy distributions and the mean absolute 
deviations from the exact Mick ipution: The optimal genera- 


tion procedure based on the mean value, m, is as follows: 


Method Mean (m) 
1. exact Poisson distribution if 0 < m < 20 
ee SU UaLe OCU transtormation -1f 20 < m < 1LO0 
Seeonnal approximation if om > 100% 
ties eu PomEOoOur DramsloOrmation should net be adopted be- 


cause it is far less accurate compared to the square root 


aly 





transformation and the normal approximation. The following 
tables and figures show the comparison statistics of the 
three approximations versus the exact Poisson distribution. 
Table I assesses the accuracy of the normal, square root, 
and cube root approximations to the exact Poisson distribu- 
Wlonmwac Selected points. Figure 1 ijimstrates graphs of the 
empirical cumulative distribution functions of the three 
approximations to that of an exact Poisson Meet iburion with 
mean 10. Figure 2 shows graphs of the mean absolute devia- 
tions of the three approximations from the exact Poisson 
CMa tivesdseterIbuUbLTOn Ilncti1om for various mean values. 


The mean absolute deviation is defined as: 


i ee : Gp = Bey? 
= a al! 


where 

Ps Lene Cb of the approximating distribution; | 

P. is the CDF of the exact Poisson distribution; and 

k is the sample size. 

Table II is a summary of a comparison of the empirical 
distributions produced by the three approximation methods To 
poicmexach FOLrsson distribution by means’ of the) one sample 
PoUmnoeorev—omnirniovy Leste. Phe null hypothesis is*®*that the 
approximated distributions are Poisson against the alterna- 
tive that they areynov Poisson. 


tasty. Wap le Itieand FPigure 3 analyze the siensitivity of 


GnessS@Uuates rooustransiormation to The constant C. in the 
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derivation of the square root transformation, C was chosen 

as 3/8. Different constants were used in order to i Neb he 
most robust constant bo Use, 1.c., LNnescOonstant which yielded 
the smallest mean absolute deviation from the exact Poisson. 
The value C = 13/18 was found to be the most robust. Note 
that this value of C was used in the approximation when 


making the comparisons with the other methods. 
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Table I. Accuracy of the Normal, Square Root, and Cube Root 
Approximations to the Exact Poisson Distribution. 








PCZ om) < n) 
square Root 3 
Normal | Transformation Cube Root 
Exact | Approx C=13/18 Transformation! 
a 0.08422 | 0.07301 Ome 7 719 0.10491 
BS 3 OF LOsvaime sl 19 39 COP IMO5 1 Omire 49 1 
6 Cpl 6 22am OReGO 36 0.14614 Om 16667 





10. |0.04861 [0.04485 0.0470 0.09045 
14 0.10244 | 0.09937 0.10288 0.05512 
18 0.07062 | 0.07622 0.07073 1 0.06348 : 
7 0.03874 | 0.03633 0.03746 | 0.08819 , 
19 0.08884 | 0.08683 0.08891 0.08726 
22 0.07692 | 0.08050 Omonieyme 0.05755 
36 0.05394 [0.05161 0.05409 0.02552 
39 0.06296 | 0.06223 0.06307 0.06336 | 
42 0.05850 | 0.05995 | 0.05841 0.05675 
1 
59 0.05019 | 0.04999 0.05020 im OOO | 
| sz 0.04007 [0.04402 0.04002 0.03289 | 
87 0.03686 | 0.03649 0.03683 0.03801 | 
89 0.03165 | 0.03243 0.03149 i 0.03991 
97 KOmosaL esp rOro3s1 82 0.03098 0.0435 
oe . 
oh 


100 
vam s09 7 | ro 39 00 0.03994 0.03605 
: 0.03627 
120 





Owe 50 
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0.014 
OiGelas 


Ceol 2 
oe Ol 
0.010 | 
2009 
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0.0064 


0.005 


mean absolute deviation 


0.004 


0.003] 
OOO 





0.001] 
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Square Root Transformation Constant. 
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IV. THE METHODS OF DIETER AND AHRENS 


During the writing of this thesis, a pre-publication 
chapter from a book by U. Dieter and J. H. Ahrens was re- 
ceived. This chapter deals with the generation of non-uni- 
form variates and contains a section on the generation of 
Poisson variates. Wie of the methods which they described 
will be outlined here. They have not been programmed or 


besved againgt #¥he composite method of this thesis. 


A. UNIT MEAN METHODS 

It is well-known that the sum of independent Poisson 
deviates with mean Uisl2.5s is Poisson distributed with mean 
HitHot... Hence the following algorithm may be considered. 

1. split the mean m of Bae Potissonvaistribucion into 
m'<+[m] and f+m-m'<1l. (Hence [m] means take the integer part 
of m and f is then the fractional part.) 

2. Take m' samples from the Poisson distribution of 
mean 1 using some suitable method. Let s be the sum of these 
samples. | 

3. Take a further sample n from the Poisson (1) distri- 
bution. If n= 0 deliver kes. Otherwise, generate U; ,U2...U, 
and count the number, j, of U, for wha cla U; Seo a eece Iver 
ks] 4 
Pic evaMicLy Ol Step 4 15 proved as follows: the probabilivy 


of observing j samples U; iesouchnoamecmmcquai torfe(once n= J 


samples greater than f) is 
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. f' (1-£) 7°. 


Since n follows a Poisson distribution with mean 1 this 











yields : 
pee bE Se (2) oS ery 2 Ee edgy td 
’ 1 ! — 4 ’ 
j an igh j | ae J'(m-j)! 
eta gue?) a eciefa Im emt 
ji (n-J5)! i! ji 


n-j=0 


The unit-mean method becomes a suitable EMbeqen oa nial fee 
the case when the exact method should not be used for gen- 
erating Poisson random variates, that is, when the mean is 
large (m > 200). This method will then become slightly faster 
then the exact method, but still needs to call too many uni- 
form deviates since it is basically the same as the exact method 
except for splitting the mean, compared to the square root 
transformation approximation method, the unit mean method may 


be uneconomical. 


B. THE CENTER TRIANGLE METHOD 

Tims =method ws suitable in the case of variable means m 
of the Poisson distribution. The mean is split into m' = [m], 
MiceeMmLeser part ofgms> and the remaining fraction Ost om-m'<1. 
Wres algorithm usdse a mixture of two dustributdions. Thesfirst 
distribution has an easy sampling method and is employed most 


Si eviceyricm~ lo Tce bDased On the following isosceles triangle. 


ef 





b 
32 


no 


g(x) = [x-m'| if m'-s < x < m'ts 


(B.1) 
0 elsewhere. 


The constant s is chosen as 
(B.2) Ss = 2,.2160358671665 vm' - 0.78125. 


The area under g(x) is b. If b is small enough then the 


inequalities 


+1 
(B.3) te veh g(x)dx < S (eS Onlie5) 


= ° 
il 
may hold for all i. The maximum b for which this is true was 
found numerically for all 4 < m < 200. The choice of the 
constant a = 2.2160358671665 in (B.2) is well motivated. The 
vertices of the maximum isosceles triangle that can be placed 
wholly inside the standard normal curve are situated at C0 s0n 


ond scam y 27) wee Asem’ Foeseto infinity, b approaches the 


max 
area of this triangle which is .88407040229876. The correc- 
tion -0.78125 = -25/32 in (B.2) improves the value of De 
for small m. [7, page 11-15, table.] Here the variable mean 
mis split into m = m'tf as in the unit mean method. If 

m < 8 the exact method is applied. If m > 9, the triangular 

method is used with the triangular probabilities te = 27/32. 

Thise@cenver triangle method needs a long)program including 


several tables which requires a large memory space and suf- 


fers@irom the’ roundoff errors' of b. 
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This method proves a very efficient algorithm based on 
a simple idea. However, it is somewhat messy to describe 
LieGmpeOnTranmands has lito les mathematical appeal compared to 
the square root transformation method. This method also 
needs a long execution time as does the exact method. For 
example, to compute one Poisson variate from a distribution 
with mean m = 100, the program would require LB66 Usml / . 


page 11-19]. 


C. THE GAMMA METHOD 

In order to obtain a sample k from the Poisson distribu- 
tion with variable mean m, select a positive integer or 
(typically n is a little smaller than m). Then, take a sam- 
ple x from a Gamma distribution with parameter n. 

Case (i) tite x > m, return a sample k from the binomial 
distribution with parameters n-l, m/x. 

Case (2) if x < m, take a sample j from the Poisson dis- 
tribution of mean m-x and return k<«ntj. 

The sample x simulates the n-th event (arrival) ina 
Poisson process of rate 1. If x > m, then there are n-l 
arrivals in the interval (0,x), and each of these has a 
CLOUD PeninveOmn, x sor being belowsm (Case h))- If x < m, 
then the n simulated arrivals are all before m and the sample 
j indicates the additional events between x and m (Case (2)). 

A Lerma eproot wot the procedure runs as follows. In the 


first case one has m < x < » and 
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-x n-l 
pe = ( n-1) (yx ete = * 


k xX X - '(n) ax 


Mm 


is the probability of obtaining k from the binomial distribu- 
tion (n-1, m/x) summed over the Gamma (n)-distributed values 


Ol yeabovye m. The expression transforms into 


p = (n-1)! m*(x—-my"71-k e7X .,n-l ae 
k a= 1.2 )2 xno n-1)! 


k co 
m Dele x 
~ ki(n-1-k)! J on) © conga 
Mm 


CO 


el 
m oe n-l-k .-t ., 
ef et 


0 


where t x —-m. Notice that r(n-k) is (n-1-k)!. Hence 


m* em 
a mee 


as required. In the second case one has 0 < x < m, and 


is the probability of obtaining j = k-n from the Poisson 
distribution of mean m-x summed over.the Gamma (n)-distri- 
buted values of x between 0 andm. This expression trans-— 


fornis. ne 
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—In 
O 


Introducing t = x/m yields 


nyse 


Be ae , nal en 


@) 


The integral is a Beta function 8(n,k-—ntl1) 


Pp = em mK T(n) l(k=n+t1) _ eM mK 
a ee Ci lwlelc=n))! T(kt1) k! 
ass betrore. The following algorithm is considered. 


ieerina Circe Ze k-O , wom 

2. If w > cC (C = 24 was used as the mean cut-off point) 
So uo. 6 

3. (Start Case (2)). Set p+1 and calculate bee” 

4, Generate U and set p+pU. If p<b deliver k 

5. Increase k+kt+1 and go to 4 

6. Set n«[dw] where d = 7/8. Take a sample x from the 
Gamma-(n) distribution. If x>w go to 8. 

{. Set k ktw, w w-x and go to 2 

8. (Start Case (1)). Set p<«w/x 

9. Generate U. If U<p inerease k<«ktl1 

LO. Set n«n-1. If n>1l go to 9 

iiewe livers ke 

The performance of the algorithm depends on the cut-off 
DOlnvecmeimabeopmes Step 3. otem 4 and Step 5 are exactly the 


Same as the steps of exact methods. This gamma method has a 


cel 








complex sub algorithm for the binomial distribution in Steps 
8-10 and requires memory space for the gamma distribution. 
This method may be efficient in cases with extremely large 
means, .say m = 1000, but this sneteincnd requires a rather long 
execution time as does the exact method. For example, to 


generate a sample from a Poisson distribution with mean 100 | 


takes 1175 us, see [7, page 11-23]. 


Ne 
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